home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / integration / qc25c.c < prev    next >
Encoding:
C/C++ Source or Header  |  2000-05-05  |  3.0 KB  |  133 lines
  1. /* integration/qc25c.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. struct fn_cauchy_params
  21. {
  22.   gsl_function *function;
  23.   double singularity;
  24. };
  25.  
  26. static double fn_cauchy (double t, void *params);
  27.  
  28. static void compute_moments (double cc, double *moment);
  29.  
  30. static void
  31. qc25c (gsl_function * f, double a, double b, double c, 
  32.        double *result, double *abserr, int *err_reliable);
  33.  
  34. static void
  35. qc25c (gsl_function * f, double a, double b, double c, 
  36.        double *result, double *abserr, int *err_reliable)
  37. {
  38.   double cc = (2 * c - b - a) / (b - a);
  39.  
  40.   if (fabs (cc) > 1.1)
  41.     {
  42.       double resabs, resasc;
  43.  
  44.       gsl_function weighted_function;
  45.       struct fn_cauchy_params fn_params;
  46.  
  47.       fn_params.function = f;
  48.       fn_params.singularity = c;
  49.  
  50.       weighted_function.function = &fn_cauchy;
  51.       weighted_function.params = &fn_params;
  52.  
  53.       gsl_integration_qk15 (&weighted_function, a, b, result, abserr,
  54.                 &resabs, &resasc);
  55.       
  56.       if (*abserr == resasc)
  57.     {
  58.       *err_reliable = 0;
  59.     }
  60.       else 
  61.     {
  62.       *err_reliable = 1;
  63.     }
  64.  
  65.       return;
  66.     }
  67.   else
  68.     {
  69.       double cheb12[13], cheb24[25], moment[25];
  70.       double res12 = 0, res24 = 0;
  71.       size_t i;
  72.       gsl_integration_qcheb (f, a, b, cheb12, cheb24);
  73.       compute_moments (cc, moment);
  74.  
  75.       for (i = 0; i < 13; i++)
  76.     {
  77.       res12 += cheb12[i] * moment[i];
  78.     }
  79.  
  80.       for (i = 0; i < 25; i++)
  81.     {
  82.       res24 += cheb24[i] * moment[i];
  83.     }
  84.  
  85.       *result = res24;
  86.       *abserr = fabs(res24 - res12) ;
  87.       *err_reliable = 0;
  88.  
  89.       return;
  90.     }
  91. }
  92.  
  93. static double
  94. fn_cauchy (double x, void *params)
  95. {
  96.   struct fn_cauchy_params *p = (struct fn_cauchy_params *) params;
  97.   gsl_function *f = p->function;
  98.   double c = p->singularity;
  99.   return GSL_FN_EVAL (f, x) / (x - c);
  100. }
  101.  
  102. static void
  103. compute_moments (double cc, double *moment)
  104. {
  105.   size_t k;
  106.  
  107.   double a0 = log (fabs ((1.0 - cc) / (1.0 + cc)));
  108.   double a1 = 2 + a0 * cc;
  109.  
  110.   moment[0] = a0;
  111.   moment[1] = a1;
  112.  
  113.   for (k = 2; k < 25; k++)
  114.     {
  115.       double a2;
  116.  
  117.       if ((k % 2) == 0)
  118.     {
  119.       a2 = 2.0 * cc * a1 - a0;
  120.     }
  121.       else
  122.     {
  123.       const double km1 = k - 1.0;
  124.       a2 = 2.0 * cc * a1 - a0 - 4.0 / (km1 * km1 - 1.0);
  125.     }
  126.  
  127.       moment[k] = a2;
  128.  
  129.       a0 = a1;
  130.       a1 = a2;
  131.     }
  132. }
  133.